Escape of a particle from an intracellular region: A particle is diffusing according to Brownian motion in a 2D circular domain of radius $r$. The particle can only exit the domain from an opening defined by a segment of a circle of angular width $\theta$.
We investigate the mean exit time of the domain as a function of the angular width $\theta$.
Brownian Motion: The particle's position follows a Brownian motion, described by the stochastic differential equation: $ d\mathbf{X}(t) = \sqrt{2D} \, d\mathbf{W}(t) $ where $ \mathbf{X}(t) $ is the position of the particle at time $ t , D $ is the diffusion coefficient, and $\mathbf{W}(t) $ is a Wiener process.
Domain: The particle moves within a circular domain of radius $ r $. The boundary conditions ensure that the particle reflects off the circle's boundary, except for an opening of angular width $ \theta $, where the particle can escape.
Exit Time: The mean exit time $ T(\theta) $ is the expected time for the particle to leave the domain through the opening.
Simulation Setup:
Visualization:
Mean Exit Time Analysis:
Histogram of Exit Times:
Density Plots:
This comprehensive study helps us understand how the angular width of the opening influences the mean exit time of particles undergoing Brownian motion within a circular domain.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
import os
import sys
import datetime
from scipy.stats import linregress
eps = 1e-6
In this section, we set up the basic simulation environment for our particle diffusion model. We fix the radius of the circle to be 1 and the diffusion coefficient to be 1. The following parameters are used for the simulation:
For the initial experiments, we set the angular width $ \theta $ to be $ \pi/4 $ to observe the general behavior of the particles.
The Experiment class defines the parameters and initial conditions for the particle diffusion simulation.
r: Radius of the circle in which the particles are diffusing. Default value is 1.0.theta: Angular width of the exit (missing arc) in the circle. Default value is π/4.D: Diffusion constant representing the rate at which particles diffuse. Default value is 1.0.dt: Time step for the simulation. Default value is 0.01.num_particles: Number of particles to simulate. Default value is 10,000.x_start: Starting x-coordinate for the particles. Default value is 0.0.y_start: Starting y-coordinate for the particles. Default value is 0.0.T_lim: Time limit for the simulation. Default value is 10,000.class Experiment:
def __init__(self, r = 1.0, theta = np.pi / 4,
D = 1.0, dt = 0.01, num_particles = 10000,
x_start = 0.0, y_start = 0.0, T_lim = 10000.):
self.r = r # radius of the circle
self.theta = theta # exit angular width
self.D = D # diffusion constant
self.dt = dt # number of particles to simulate
self.num_particles = num_particles # number of particles to simulate
self.x_start = x_start # starting x coordinate
self.y_start = y_start # starting y coordinate
self.T_lim = T_lim
E = Experiment
# Parameters
default_experiment = E(
r = 1.0,
theta = np.pi / 4,
D = 1.0,
dt = 0.01,
num_particles = 10000,
x_start = 0.0,
y_start = 0.0
)
The find_intersection function determines the intersection point between the path of a particle and the boundary of a circle. This is useful for simulating particle movement and checking if a particle has escaped the circular region.
x: Current x-coordinate of the particle.y: Current y-coordinate of the particle.x_prev: Previous x-coordinate of the particle.y_prev: Previous y-coordinate of the particle.r: Radius of the circle.(x_int, y_int): Coordinates of the intersection point if it exists. If no intersection is found within the valid range, it returns (None, None).def find_intersection(x, y, x_prev, y_prev, r):
x_prev_norm_squared = x_prev**2 + y_prev**2
x_norm_squared = x**2 + y**2
dot_product = x_prev * x + y_prev * y
a = x_prev_norm_squared + x_norm_squared - 2 * dot_product
b = 2 * dot_product - 2 * x_norm_squared
c = x_norm_squared - r**2
discriminant = b**2 - 4 * a * c
if discriminant < 0:
return None, None
lambda_1 = (-b + np.sqrt(b**2 - 4 * a * c)) / (2 * a)
lambda_2 = (-b - np.sqrt(b**2 - 4 * a * c)) / (2 * a)
if 0 < lambda_1 < 1:
x_int = x + lambda_1 * (x_prev - x)
y_int = y + lambda_1 * (y_prev - y)
elif 0 < lambda_2 < 1:
x_int = x + lambda_2 * (x_prev - x)
y_int = y + lambda_2 * (y_prev - y)
else:
return None, None
return x_int, y_int
The is_in_exit_segment function determines if the intersection point between a particle's path and the circle boundary lies within the defined exit segment (missing arc).
x: Current x-coordinate of the particle.y: Current y-coordinate of the particle.x_prev: Previous x-coordinate of the particle.y_prev: Previous y-coordinate of the particle.theta: Angular width of the exit segment (missing arc).r: Radius of the circle.bool: True if the intersection point is within the exit segment, False otherwise.def is_in_exit_segment(x, y, x_prev, y_prev, theta, r):
# find intersection point
x_int, y_int = find_intersection(x, y, x_prev, y_prev, r)
if x_int is None:
return False
# check if intersection point is in the exit segment
return -theta/2 < np.arctan2(y_int, x_int) < theta / 2
These functions handle the reflection of particles and their Brownian motion within a circular boundary with a missing arc.
The reflect function calculates the new position of a particle after it reflects off the boundary of the circle.
x: Current x-coordinate of the particle.y: Current y-coordinate of the particle.x_prev: Previous x-coordinate of the particle.y_prev: Previous y-coordinate of the particle.r: Radius of the circle.(x_reflection, y_reflection): Coordinates of the particle after reflection.def reflect(x, y, x_prev, y_prev, r):
# Calculate the coefficients of the quadratic equation
x_intersection, y_intersection = find_intersection(x, y, x_prev, y_prev, r)
if x_intersection is None:
raise ValueError("No intersection point found")
# Calculate the normal vector at the intersection point
normal_x = x_intersection / r
normal_y = y_intersection / r
# Calculate the incident vector
incident_x = x - x_intersection
incident_y = y - y_intersection
# Reflect the incident vector across the normal vector
dot_product = incident_x * normal_x + incident_y * normal_y
reflection_x = incident_x - 2 * dot_product * normal_x
reflection_y = incident_y - 2 * dot_product * normal_y
# Calculate the reflected position
x_reflection = x_intersection + reflection_x
y_reflection = y_intersection + reflection_y
return x_reflection, y_reflection
brownian_step¶This function simulates a single Brownian step for a particle within the circular region and handles its reflection if it hits the boundary.
x: Current x-coordinate of the particle.y: Current y-coordinate of the particle.experiment: An instance of the Experiment class containing simulation parameters.(x, y): New coordinates of the particle after the Brownian step.def brownian_step(x, y, experiment):
x_prev, y_prev = x, y
D = experiment.D
dt = experiment.dt
theta = experiment.theta
r = experiment.r
# Brownian motion: update position with random steps
step_size = np.sqrt(2 * D * dt)
dx = step_size * np.random.randn()
dy = step_size * np.random.randn()
x_new = x + dx
y_new = y + dy
# Check if particle is outside the circle
if x_new**2 + y_new**2 > r**2:
if is_in_exit_segment(x_new, y_new, x_prev, y_prev, theta, r):
x, y = x_new, y_new
else:
# Reflect the particle back
try:
x, y = reflect(x_new, y_new, x, y, r)
except ValueError:
# Handle the case where reflection calculation fails
x, y = x_prev, y_prev # Reset to previous position
else:
x, y = x_new, y_new
return x, y
brownian_step_escaped¶This function simulates a single Brownian step for a particle that has escaped the circular region and handles its reflection if it hits the boundary.
x: Current x-coordinate of the particle.y: Current y-coordinate of the particle.experiment: An instance of the Experiment class containing simulation parameters.(x, y): New coordinates of the particle after the Brownian step.def brownian_step_escaped(x, y, experiment):
x_prev, y_prev = x, y
D = experiment.D
dt = experiment.dt
theta = experiment.theta
r = experiment.r
# Brownian motion: update position with random steps
step_size = np.sqrt(2 * D * dt)
dx = step_size * np.random.randn()
dy = step_size * np.random.randn()
x_new = x + dx
y_new = y + dy
# Check if particle is outside the circle
if x_new**2 + y_new**2 < r**2 + eps:
if is_in_exit_segment(x_new, y_new, x, y, theta, r):
x, y = x_new, y_new
else:
# Reflect the particle back
try:
x, y = reflect(x_new, y_new, x, y, r)
except ValueError:
# Handle the case where reflection calculation fails
x, y = x_prev, y_prev # Reset to previous position
else:
x, y = x_new, y_new
return x, y
Here we plot a graph to display how our reflect function works. We can see that the particle is reflected off the boundary of the circle.
def reflect_func(experiment):
x_start = 0.5
y_start = 0.5
x_new = -0.5
y_new = 1.5
x_int, y_int = find_intersection(x_new, y_new, x_start, y_start, experiment.r)
print(x_int, y_int)
x_reflected, y_reflected = reflect(x_new, y_new, x_start, y_start, experiment.r)
fig, ax = plt.subplots(figsize=(15, 10))
# mark particle points, write label next to them
ax.plot(x_start, y_start, 'ro', label='Starting point')
ax.text(x_start, y_start, 'Starting \npoint', fontsize=12, ha='left')
ax.plot(x_new, y_new, 'ro', label='New point')
ax.text(x_new, y_new, 'New point', fontsize=12, ha='right')
ax.plot(x_int, y_int, 'ro', label='Intersection point')
ax.text(x_int, y_int, 'Intersection\n point', fontsize=12, ha='left')
ax.plot(x_reflected, y_reflected, 'ro', label='Reflected point')
ax.text(x_reflected, y_reflected, 'Reflected\n point', fontsize=12, ha='right')
# draw line between points
ax.plot([x_start, x_new], [y_start, y_new], 'b-', label='Line between starting and new point')
ax.plot([x_int, x_reflected], [y_int, y_reflected], 'g-', label='Line between intersection and reflected point')
# draw arrows in the middle of the lines
ax.annotate('', xy=(0.5*x_start + 0.5*x_int, 0.5*y_start + 0.5*y_int), xytext=(x_start, y_start), arrowprops=dict(arrowstyle='->', color='b'))
ax.annotate('', xy=(0.5*x_int + 0.5*x_new, 0.5*y_int + 0.5*y_new), xytext=(x_start, y_start), arrowprops=dict(arrowstyle='->', color='b'))
ax.annotate('', xy=(0.5*x_reflected + 0.5*x_int, 0.5*y_reflected + 0.5*y_int), xytext=(x_int, y_int), arrowprops=dict(arrowstyle='->', color='g'))
# draw tangent line at intersection point
ax.plot([x_int - y_int, x_int + y_int], [y_int + x_int, y_int - x_int], 'g--', label='Tangent line at intersection point')
theta = experiment.theta
r = experiment.r
margin = 0.1 * r
# Define the circle with a missing arc
theta_start = theta / 2. # Start angle for the arc (in radians)
theta_end = 2 * np.pi - theta / 2. # End angle for the arc (in radians)
# Create the solid arc
arc_solid = np.linspace(theta_start, theta_end, 100)
circle_x_solid = r * np.cos(arc_solid)
circle_y_solid = r * np.sin(arc_solid)
# Create the dotted arc for the missing segment
arc_dotted = np.linspace(theta_end, 2 * np.pi + theta_start, 100)
circle_x_dotted = r * np.cos(arc_dotted)
circle_y_dotted = r * np.sin(arc_dotted)
# Plot the solid arc
ax.plot(circle_x_solid, circle_y_solid, color='y')
# Plot the dotted arc
ax.plot(circle_x_dotted, circle_y_dotted, 'g--')
# Set the aspect of the plot to be equal
ax.set_aspect('equal')
# Set the limits to create a box from (-r, -r) to (r, r)
ax.set_xlim(-r - margin, 2*r + margin)
ax.set_ylim(-r - margin, 2*r + margin)
# Add grid lines for better visualization
ax.grid(True)
ax.legend()
# Display the plot
plt.show()
reflect_func(E())
0.0 1.0
The simulate_particle function simulates the movement of a single particle within a circular region with a missing arc, tracking the time it takes for the particle to escape the circle.
experiment: An instance of the Experiment class containing simulation parameters.time: The time it takes for the particle to escape the circle. If the particle does not escape within the time limit, it returns None.def simulate_particle(experiment):
# Start at the center
x, y = experiment.x_start, experiment.y_start
time = 0.0
D = experiment.D
dt = experiment.dt
theta = experiment.theta
r = experiment.r
while time < experiment.T_lim:
x, y = brownian_step(x, y, experiment)
time += dt
if x**2 + y**2 > r**2:
return time
return None
We write a function to plot the random paths of a particle in the circle.
def plot_path(experiment):
theta = experiment.theta
r = experiment.r
margin = 0.5 * r
# Create the plot
fig, ax = plt.subplots(figsize=(10, 10))
plt.title(f'Path of Particles in a Circle with a Missing Arc, starting point = {experiment.x_start, experiment.y_start}')
plt.xlabel('x')
plt.ylabel('y')
# Plot the path of the particles
for _ in range(experiment.num_particles):
x, y = experiment.x_start, experiment.y_start
x_prev, y_prev = x, y
x_path = [x]
y_path = [y]
time = 0.0
while time < experiment.T_lim:
x, y = brownian_step(x, y, experiment)
x_path.append(x)
y_path.append(y)
time += experiment.dt
if x**2 + y**2 > r**2:
break
plt.plot(x_path, y_path, label=f'Path of Particle{_+1}')
# Define the circle with a missing arc
theta_start = theta / 2. # Start angle for the arc (in radians)
theta_end = 2 * np.pi - theta / 2. # End angle for the arc (in radians)
# Create the solid arc
arc_solid = np.linspace(theta_start, theta_end, 100)
circle_x_solid = r * np.cos(arc_solid)
circle_y_solid = r * np.sin(arc_solid)
# Create the dotted arc for the missing segment
arc_dotted = np.linspace(theta_end, 2 * np.pi + theta_start, 100)
circle_x_dotted = r * np.cos(arc_dotted)
circle_y_dotted = r * np.sin(arc_dotted)
# Plot the solid arc
ax.plot(circle_x_solid, circle_y_solid, color='y', label=f'Circle with radius {r} and missing arc of {theta * 180 / np.pi:.1f}°')
# Plot the dotted arc
ax.plot(circle_x_dotted, circle_y_dotted, 'g--', label='Missing arc')
# Set the aspect of the plot to be equal
ax.set_aspect('equal')
# Set the limits to create a box from (-r, -r) to (r, r)
ax.set_xlim(-r - margin, r + margin)
ax.set_ylim(-r - margin, r + margin)
# Add grid lines for better visualization
ax.grid(True)
#highlight starting point
ax.plot(experiment.x_start, experiment.y_start, 'ro', label='Starting point')
ax.legend()
# Display the plot
plt.show()
# Plot the path of a particle
experiment = E(theta=np.pi/4, num_particles=2, D=1.0)
plot_path(experiment)
We plot the random path of the particle also starting at (0.5,0.5) to see how the particle moves.
experiment = E(theta=np.pi/4, num_particles=2, D=1.0, x_start=0.5, y_start=0.5)
plot_path(experiment)
We simulate the mean exit time for a large number of particles to understand their behavior within a circle with a missing arc.
mean_exit_times = [simulate_particle(default_experiment) for _ in range(default_experiment.num_particles)]
mean_exit_time = np.mean(mean_exit_times)
print(f"Mean exit time (discretized with Brownian motion and elastic reflection): {mean_exit_time}")
Mean exit time (discretized with Brownian motion and elastic reflection): 2.2679979999999818
# Plot histogram of exit times
plt.hist(mean_exit_times, bins=150)
plt.xlabel('Exit Time')
plt.ylabel('Frequency')
plt.title('Histogram of Exit Times (Discretized with Brownian motion and elastic reflection)')
plt.show()
Now, we alter the parameter $\theta$ to see how the mean exit time changes.
def mean_exit_time_theta(experiment):
x_start = experiment.x_start
y_start = experiment.y_start
theta = experiment.theta
num_particles = experiment.num_particles
mean_exit_times = [simulate_particle(E(theta=theta, x_start=x_start, y_start=y_start)) for _ in range(num_particles)]
mean_exit_time = np.mean(mean_exit_times)
return mean_exit_time
# run simulation with different theta values
theta_values = np.linspace(np.pi/16, 2*np.pi, 30)
start_time = time.time()
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta = theta, num_particles=10000)) for theta in theta_values]
end_time = time.time()
print(f"Time taken to run simulation with different theta values: {end_time - start_time} seconds")
# Plot mean exit time as a function of theta
print(mean_exit_times_theta_list)
print(theta_values)
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.show()
Time taken to run simulation with different theta values: 142.19080805778503 seconds [4.801980999999962, 3.276215999999968, 2.625478999999976, 2.191859999999982, 1.9087109999999865, 1.6736339999999899, 1.5094939999999921, 1.3455059999999939, 1.2263979999999954, 1.1231749999999963, 1.0183389999999974, 0.9103399999999984, 0.8408159999999987, 0.7828739999999992, 0.7101369999999995, 0.6582799999999996, 0.6040289999999998, 0.5689699999999999, 0.5246970000000001, 0.4810530000000002, 0.4505780000000001, 0.42038600000000015, 0.3868130000000002, 0.36856100000000014, 0.34138700000000016, 0.33553700000000014, 0.3184230000000001, 0.3062100000000001, 0.2978460000000001, 0.29323700000000014] [0.19634954 0.40624043 0.61613132 0.82602221 1.03591309 1.24580398 1.45569487 1.66558576 1.87547665 2.08536754 2.29525843 2.50514931 2.7150402 2.92493109 3.13482198 3.34471287 3.55460376 3.76449465 3.97438553 4.18427642 4.39416731 4.6040582 4.81394909 5.02383998 5.23373086 5.44362175 5.65351264 5.86340353 6.07329442 6.28318531]
The mean exit time seems to decay at an exponential-like rate. The asymptote is well-defined beyond $\theta > 2\pi$, indicating no boundaries. As $\theta$ approaches $0$, the mean exit time should go to $\rightarrow \infty$ since the particles never leave the circle.
Here, we not only alter the parameter $\theta$, but also the starting position of the particle. We will see how the mean exit time changes with the starting position of the particle.
# change starting position of particles
mean_exit_times_theta_list = []
start_time = time.time()
for theta in theta_values:
experiment = Experiment(theta = theta, num_particles = 10000, x_start = 0.5, y_start = 0.5)
mean_exit_times_theta_list.append(mean_exit_time_theta(experiment))
end_time = time.time()
print(f"Time taken to run simulation with different theta values and starting position = (0.5,0.5): {end_time - start_time} seconds")
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta (Starting Position: x=0.5, y=0.5)')
plt.show()
Time taken to run simulation with different theta values and starting position = (0.5,0.5): 97.54851698875427 seconds
# plot both curves
theta_values = np.linspace(np.pi/16, 2*np.pi, 30)
mean_exit_times_theta_list_moved = [mean_exit_time_theta(E(theta = theta, num_particles=10000, x_start=0.5, y_start=0.5)) for theta in theta_values]
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta=theta, num_particles=10000)) for theta in theta_values]
plt.plot(theta_values, mean_exit_times_theta_list, marker='o', label='Starting Position: x=0, y=0')
plt.plot(theta_values, mean_exit_times_theta_list_moved, marker='o', label='Starting Position: x=0.5, y=0.5')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.legend()
plt.show()
Intuitively, we expect the mean exit time for the particle starting at (0.5,0.5) will be lower because it is closer to the opening.
Empirically, we observe a decay, and it is interesting to see how the rate depends on the size of the exit segment. We tried fitting different functions, and the best fit was an inverse linear function.
We try to fit this function to the mean exit time function \begin{equation} T(\theta) = \beta e^{\alpha / \theta} \end{equation}
We chose this one because of the exponential decay and clearly defined horizontal asymptote like the one we observe in our graphs.
def fit_min_exit_time_to_exp(theta_vals, mean_exit_time_vals):
my_theta_vals = np.array(theta_vals)
my_mean_exit_time_vals = np.array(mean_exit_time_vals)
x = 1 / my_theta_vals
y = np.log(my_mean_exit_time_vals)
# Step 3: Use linregress to find the slope and intercept
slope, intercept, r_value, p_value, std_err = linregress(x, y)
alpha = slope
beta = np.exp(intercept)
return alpha, beta
# run simulation with different theta values
theta_values = np.linspace(np.pi/16, 2*np.pi, 30)
start_time = time.time()
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta = theta, num_particles=10000)) for theta in theta_values]
alpha, beta = fit_min_exit_time_to_exp(theta_values, mean_exit_times_theta_list)
print(alpha, beta)
x = np.linspace(np.pi/16, 2*np.pi, 100)
plt.plot(x, beta * np.exp(beta / x),
'r-', lw=5, alpha=0.6, label='levy pdf')
end_time = time.time()
print(f"Time taken to run simulation with different theta values: {end_time - start_time} seconds")
# Plot mean exit time as a function of theta
print(mean_exit_times_theta_list)
print(theta_values)
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.show()
0.6083530268701564 0.5192810677166473 Time taken to run simulation with different theta values: 49.34588432312012 seconds [4.652655999999959, 3.244538999999967, 2.5632949999999766, 2.175579999999983, 1.881732999999987, 1.6848109999999898, 1.4841629999999923, 1.3448729999999942, 1.227317999999995, 1.1157899999999965, 1.0043569999999975, 0.9351009999999981, 0.8392119999999986, 0.7607329999999991, 0.7069919999999995, 0.6456379999999997, 0.6023369999999999, 0.5537040000000001, 0.5071490000000001, 0.47117300000000023, 0.44061900000000026, 0.41995000000000016, 0.39148200000000016, 0.3680220000000001, 0.3490230000000002, 0.33379400000000015, 0.31866900000000004, 0.30339300000000013, 0.3007040000000001, 0.2958450000000001] [0.19634954 0.40624043 0.61613132 0.82602221 1.03591309 1.24580398 1.45569487 1.66558576 1.87547665 2.08536754 2.29525843 2.50514931 2.7150402 2.92493109 3.13482198 3.34471287 3.55460376 3.76449465 3.97438553 4.18427642 4.39416731 4.6040582 4.81394909 5.02383998 5.23373086 5.44362175 5.65351264 5.86340353 6.07329442 6.28318531]
Although it has a horizontal asymptote, the inverse linear function has a smaller decay rate. It turns out to be a pretty good fit for the observed data. \begin{equation} T(\theta) = \beta + \frac{\alpha}{\theta} \end{equation}
from scipy.optimize import curve_fit
def fit_min_exit_time_to_inv_linear(theta_vals, mean_exit_time_vals):
my_theta_vals = np.array(theta_vals)
my_mean_exit_time_vals = np.array(mean_exit_time_vals)
x = 1 / my_theta_vals
y = my_mean_exit_time_vals
# Step 3: Use linregress to find the slope and intercept
slope, intercept, r_value, p_value, std_err = linregress(x, y)
alpha = slope
beta = intercept
return alpha, beta
# run simulation with different theta values
theta_values = np.linspace(np.pi/16, 2*np.pi, 30)
start_time = time.time()
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta = theta, num_particles=10000)) for theta in theta_values]
alpha, beta = fit_min_exit_time_to_inv_linear(theta_values, mean_exit_times_theta_list)
print(alpha, beta)
x = np.linspace(np.pi/16, 2*np.pi, 100)
plt.plot(x, beta + alpha / x,
'r-', lw=5, alpha=0.6, label='levy pdf')
end_time = time.time()
print(f"Time taken to run simulation with different theta values: {end_time - start_time} seconds")
# Plot mean exit time as a function of theta
print(mean_exit_times_theta_list)
print(theta_values)
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.show()
0.9603724566633827 0.43373437871169873 Time taken to run simulation with different theta values: 49.462323904037476 seconds [4.630552999999961, 3.1502539999999692, 2.5508549999999772, 2.1572829999999827, 1.9160469999999867, 1.6769069999999902, 1.4898019999999923, 1.3394319999999942, 1.2277189999999951, 1.0880839999999965, 0.9962579999999974, 0.9150159999999982, 0.8438809999999987, 0.7687909999999992, 0.7275759999999994, 0.6497929999999996, 0.6044989999999999, 0.552826, 0.52054, 0.48250200000000015, 0.4421420000000002, 0.4059690000000002, 0.39083500000000015, 0.36682000000000015, 0.3441730000000001, 0.33121000000000017, 0.32054500000000014, 0.3051520000000002, 0.30050100000000013, 0.29565300000000005] [0.19634954 0.40624043 0.61613132 0.82602221 1.03591309 1.24580398 1.45569487 1.66558576 1.87547665 2.08536754 2.29525843 2.50514931 2.7150402 2.92493109 3.13482198 3.34471287 3.55460376 3.76449465 3.97438553 4.18427642 4.39416731 4.6040582 4.81394909 5.02383998 5.23373086 5.44362175 5.65351264 5.86340353 6.07329442 6.28318531]
Overall, the inverse linear function offers a more accurate representation of the mean exit time as a function of $\theta$, balancing both the decay rate and the asymptotic behavior effectively.
curve_fit Method for Better Fitting¶Let's use the curve_fit method instead of manual fitting to achieve a more accurate fit for the decay function. This approach leverages the optimization capabilities of curve_fit to fine-tune the parameters for the best fit.
def fit_min_exit_time_to_inv_linear_with_curvefit(theta_vals, mean_exit_time_vals):
def formula(x, alpha, beta):
return beta + alpha / x
params, covariance = curve_fit(formula, theta_vals, mean_exit_time_vals)
alpha, beta = params
return alpha, beta
# run simulation with different theta values
theta_values = np.linspace(np.pi/16, 2*np.pi, 30)
start_time = time.time()
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta = theta, num_particles=10000)) for theta in theta_values]
alpha, beta = fit_min_exit_time_to_inv_linear_with_curvefit(theta_values, mean_exit_times_theta_list)
print(alpha, beta)
x = np.linspace(np.pi/16, 2*np.pi, 100)
plt.plot(x, beta + alpha / x,
'r-', lw=5, alpha=0.6, label='levy pdf')
end_time = time.time()
print(f"Time taken to run simulation with different theta values: {end_time - start_time} seconds")
# Plot mean exit time as a function of theta
print(mean_exit_times_theta_list)
print(theta_values)
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.show()
0.9578588725370406 0.4343561065718333 Time taken to run simulation with different theta values: 46.6053671836853 seconds [4.609769999999962, 3.159815999999968, 2.5647269999999773, 2.159324999999983, 1.8632929999999872, 1.6824189999999901, 1.5045399999999922, 1.358819999999994, 1.2111269999999956, 1.1176969999999964, 0.9915279999999977, 0.9014269999999985, 0.8387399999999987, 0.767757999999999, 0.7138259999999995, 0.6496969999999996, 0.6089059999999997, 0.559756, 0.5181720000000001, 0.4746750000000002, 0.44567300000000015, 0.41586500000000015, 0.38876100000000013, 0.36105700000000013, 0.34853900000000015, 0.32718600000000014, 0.3141080000000001, 0.30303800000000014, 0.30110800000000015, 0.29976400000000014] [0.19634954 0.40624043 0.61613132 0.82602221 1.03591309 1.24580398 1.45569487 1.66558576 1.87547665 2.08536754 2.29525843 2.50514931 2.7150402 2.92493109 3.13482198 3.34471287 3.55460376 3.76449465 3.97438553 4.18427642 4.39416731 4.6040582 4.81394909 5.02383998 5.23373086 5.44362175 5.65351264 5.86340353 6.07329442 6.28318531]
This approach aims to reduce overfitting and improve the generalization of the fitted function to accurately describe the mean exit time as a function of $\theta$.
# run simulation with different theta values
L = 8 *np.pi
theta_values = np.linspace(np.pi/16, L, 50)
start_time = time.time()
mean_exit_times_theta_list = [mean_exit_time_theta(E(theta = theta, num_particles=10000)) for theta in theta_values]
alpha, beta = fit_min_exit_time_to_inv_linear_with_curvefit(theta_values, mean_exit_times_theta_list)
print(alpha, beta)
x = np.linspace(np.pi/16, L, 100)
plt.plot(x, beta + alpha / x,
'r-', lw=5, alpha=0.6, label='levy pdf')
end_time = time.time()
print(f"Time taken to run simulation with different theta values: {end_time - start_time} seconds")
# Plot mean exit time as a function of theta
print(mean_exit_times_theta_list)
print(theta_values)
plt.plot(theta_values, mean_exit_times_theta_list, marker='o')
plt.xlabel('Theta')
plt.ylabel('Mean Exit Time')
plt.title('Mean Exit Time as a Function of Theta')
plt.show()
0.9393020031843846 0.2755077730370878 Time taken to run simulation with different theta values: 36.72853183746338 seconds [4.680448999999961, 2.35597799999998, 1.7114829999999892, 1.3117439999999945, 1.0337549999999973, 0.8473549999999989, 0.6761969999999997, 0.565031, 0.47027400000000014, 0.3995010000000002, 0.3416500000000001, 0.30870600000000015, 0.2942590000000001, 0.2953540000000001, 0.2941470000000001, 0.2964810000000001, 0.2926560000000001, 0.29368200000000017, 0.29418100000000014, 0.2932300000000001, 0.2947230000000001, 0.29391100000000003, 0.29388700000000006, 0.29407300000000014, 0.2957760000000001, 0.29615900000000006, 0.29306900000000013, 0.29628000000000004, 0.2981710000000001, 0.2952000000000002, 0.29435300000000014, 0.29304400000000014, 0.29749600000000015, 0.29730600000000007, 0.2950570000000002, 0.2943300000000001, 0.29487600000000014, 0.29607900000000015, 0.2952470000000001, 0.29657400000000017, 0.3025850000000001, 0.29606800000000016, 0.2947410000000001, 0.29416400000000015, 0.2961660000000001, 0.2899330000000001, 0.29568100000000014, 0.29501700000000014, 0.2961630000000001, 0.2936220000000001] [ 0.19634954 0.70525549 1.21416145 1.7230674 2.23197335 2.7408793 3.24978526 3.75869121 4.26759716 4.77650312 5.28540907 5.79431502 6.30322097 6.81212693 7.32103288 7.82993883 8.33884479 8.84775074 9.35665669 9.86556264 10.3744686 10.88337455 11.3922805 11.90118646 12.41009241 12.91899836 13.42790431 13.93681027 14.44571622 14.95462217 15.46352813 15.97243408 16.48134003 16.99024598 17.49915194 18.00805789 18.51696384 19.02586979 19.53477575 20.0436817 20.55258765 21.06149361 21.57039956 22.07930551 22.58821146 23.09711742 23.60602337 24.11492932 24.62383528 25.13274123]
We see that using a wider ranges of $\theta$-s give us a quite statisfying approximation for our data. The approximation is not perfect and it is probable that the underlying function is more complex, however we stop our research here.
The histogram represents the distribution of exit times for particles with a fixed $\theta$. It visualizes how frequently particles escape the circular region at different times.
The Levy distribution is a stable distribution, and one of the few that is stable under addition. This makes it suitable for modeling the time between events in Brownian motion, as it accurately represents the heavy-tailed nature of the distribution of exit times.
mean_exit_times_tmp = [simulate_particle(E(r = 1., T_lim = 20.)) for _ in range(10000)]
mean_exit_times_tmp = [val for val in mean_exit_times_tmp if val is not None]
mean_exit_time_tmp = np.mean(mean_exit_times)
print(f"Mean exit time (discretized with Brownian motion and elastic reflection): {mean_exit_time}")
Mean exit time (discretized with Brownian motion and elastic reflection): 2.2279219999999818
# Plot histogram of exit times
plt.hist(mean_exit_times_tmp, bins=30)
plt.xlabel('Exit Time')
plt.ylabel('Frequency')
plt.title('Histogram of Exit Times (Discretized with Brownian motion and elastic reflection)')
plt.show()
from matplotlib.pylab import norm
from scipy.stats import levy
loc, scale = levy.fit(mean_exit_times_tmp)
print(scale, loc)
x = np.linspace(0, 10, 100)
plt.xlim(0, 10)
plt.plot(x, levy.pdf(x, loc = loc, scale = scale),
'r-', lw=5, alpha=0.6, label='levy pdf')
plt.hist(mean_exit_times_tmp, bins=150, density = True)
plt.legend()
plt.show()
0.7189944278186888 -0.009977055599261802
In conclusion, the Levy distribution's heavy tail and flexibility make it an excellent fit for modeling the exit times of particles undergoing Brownian motion, capturing both the initial rise and the long tail observed in the empirical data.
These density plots provide a visual understanding of the diffusion dynamics of particles, highlighting how their distribution evolves over time.
# In parallel simulation basically and scatter plot at each step.
from re import I
def make_snapshot(experiment, coords, escaped_coords, mode = 'scatter'):
# Parameters
r = experiment.r # Radius of the circle
theta = experiment.r # Angle of the missing arc in radians (e.g., 45 degrees)
margin = 0.5
# Get x and y
x = np.array([_[0] for _ in coords])
y = np.array([_[1] for _ in coords])
x_esc = np.array([_[0] for _ in escaped_coords])
y_esc = np.array([_[1] for _ in escaped_coords])
# Create the plot
fig, ax = plt.subplots()
if mode == 'scatter':
# Scatter plot of the points
ax.scatter(x, y, label='Data points')
ax.scatter(x_esc, y_esc, label='Escaped data points', color = 'red')
else:
x_full = np.concatenate((x, x_esc))
y_full = np.concatenate((y, y_esc))
# Create the heatmap
heatmap, xedges, yedges = np.histogram2d(x_full, y_full, bins=50, range=[[-r - margin, r + margin], [-r - margin, r + margin]])
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
# Plot the heatmap
cax = ax.imshow(heatmap.T, extent=extent, origin='lower', cmap='viridis', aspect='auto')
# Add a colorbar
cbar = fig.colorbar(cax, ax=ax)
cbar.set_label('Frequency')
# Define the circle with a missing arc
theta_start = theta / 2. # Start angle for the arc (in radians)
theta_end = 2 * np.pi - theta / 2. # End angle for the arc (in radians)
# Create the solid arc
arc_solid = np.linspace(theta_start, theta_end, 100)
circle_x_solid = r * np.cos(arc_solid)
circle_y_solid = r * np.sin(arc_solid)
# Create the dotted arc for the missing segment
arc_dotted = np.linspace(theta_end, 2 * np.pi + theta_start, 100)
circle_x_dotted = r * np.cos(arc_dotted)
circle_y_dotted = r * np.sin(arc_dotted)
# Plot the solid arc
ax.plot(circle_x_solid, circle_y_solid, color='y', label=f'Circle with radius {r} and missing arc of {theta * 180 / np.pi:.1f}°')
# Plot the dotted arc
ax.plot(circle_x_dotted, circle_y_dotted, 'g--', label='Missing arc')
# Set the aspect of the plot to be equal
ax.set_aspect('equal')
# Set the limits to create a box from (-r, -r) to (r, r)
ax.set_xlim(-r - margin, r + margin)
ax.set_ylim(-r - margin, r + margin)
# Add grid lines for better visualization
ax.grid(True)
# Add a legend
ax.legend()
# Display the plot
plt.show()
def simulate_particle_with_snapshots(experiment, snapshot_steps, mode = 'scatter'):
# Start at the center
coords = [(experiment.x_start, experiment.y_start)] * experiment.num_particles
escaped_coords = []
time = 0.0
step = 0
max_step = np.max(snapshot_steps)
D = experiment.D
dt = experiment.dt
theta = experiment.theta
r = experiment.r
while step < max_step:
new_coords = []
new_escaped_coords = []
for (i, (x, y)) in enumerate(escaped_coords):
x, y = brownian_step_escaped(x, y, experiment)
escaped_coords[i] = (x, y)
for (x, y) in coords:
x, y = brownian_step(x, y, experiment)
if x**2 + y**2 <= r**2:
new_coords.append((x, y))
else:
escaped_coords.append((x, y))
time += dt
step += 1
coords = new_coords
if step in snapshot_steps:
print(f'Snapshot at time t = {round(time, 3)}:')
make_snapshot(experiment, coords, escaped_coords, mode)
return None
simulate_particle_with_snapshots(E(num_particles=100), range(1, 100, 5))
Snapshot at time t = 0.01:
Snapshot at time t = 0.06:
Snapshot at time t = 0.11:
Snapshot at time t = 0.16:
Snapshot at time t = 0.21:
Snapshot at time t = 0.26:
Snapshot at time t = 0.31:
Snapshot at time t = 0.36:
Snapshot at time t = 0.41:
Snapshot at time t = 0.46:
Snapshot at time t = 0.51:
Snapshot at time t = 0.56:
Snapshot at time t = 0.61:
Snapshot at time t = 0.66:
Snapshot at time t = 0.71:
Snapshot at time t = 0.76:
Snapshot at time t = 0.81:
Snapshot at time t = 0.86:
Snapshot at time t = 0.91:
Snapshot at time t = 0.96:
simulate_particle_with_snapshots(E(num_particles=10000), range(1, 100, 5), mode = 'heatmap')
Snapshot at time t = 0.01:
Snapshot at time t = 0.06:
Snapshot at time t = 0.11:
Snapshot at time t = 0.16:
Snapshot at time t = 0.21:
Snapshot at time t = 0.26:
Snapshot at time t = 0.31:
Snapshot at time t = 0.36:
Snapshot at time t = 0.41:
Snapshot at time t = 0.46:
Snapshot at time t = 0.51:
Snapshot at time t = 0.56:
Snapshot at time t = 0.61:
Snapshot at time t = 0.66:
Snapshot at time t = 0.71:
Snapshot at time t = 0.76:
Snapshot at time t = 0.81:
Snapshot at time t = 0.86:
Snapshot at time t = 0.91:
Snapshot at time t = 0.96:
From a position (0.5, 0.5)
simulate_particle_with_snapshots(E(num_particles=10000, x_start = 0.5, y_start = 0.5), range(1, 100, 5), mode = 'heatmap')
Snapshot at time t = 0.01:
Snapshot at time t = 0.06:
Snapshot at time t = 0.11:
Snapshot at time t = 0.16:
Snapshot at time t = 0.21:
Snapshot at time t = 0.26:
Snapshot at time t = 0.31:
Snapshot at time t = 0.36:
Snapshot at time t = 0.41:
Snapshot at time t = 0.46:
Snapshot at time t = 0.51:
Snapshot at time t = 0.56:
Snapshot at time t = 0.61:
Snapshot at time t = 0.66:
Snapshot at time t = 0.71:
Snapshot at time t = 0.76:
Snapshot at time t = 0.81:
Snapshot at time t = 0.86:
Snapshot at time t = 0.91:
Snapshot at time t = 0.96:
We appreciate your time and hope you found the insights valuable.
For more information or further inquiries, feel free to reach out.
Happy Analyzing! 📊